C ***********************************************************************
C *            Subroutine LIBEMA by Stephen Kirkup                      *
C ***********************************************************************
C
C Website : www.boundary-element-method.com (LAPLACE/BEMLAP)
C
C This subroutine computes the solution to the axisymmetric Laplace
C equation
C                  __ 2           
C                  \/    {\phi}    =  0   
C
C in the domain interior to a closed axisymmetric boundary.
C
C The boundary (S) is defined (approximated) by a set of conical 
C panels. The domain of the equation is within the boundary.
C
C The boundary condition may be Dirichlet, Robin or Neumann. It is also
C axisymmetric and it is assumed to have the following general form
C
C           {\alpha}(q) {\phi}(q) + {\beta}(q) v(q) = f(q)
C    
C where {\phi}(q) is the solution at the point q on S, v(q) is the 
C derivative of {\phi} with respect to the outward normal to S at q and
C {\alpha}, {\beta} and f are real-valued functions defined on S. The
C functions {\alpha} and {\beta} must be specified to define the nature
C of the boundary condition. Important examples are {\alpha}=1, 
C {\beta}=0 which is equivalent to a Dirichlet boundary condition and 
C {\alpha}=0, {\beta}=1 which is equivalent to a Neumann boundary 
C condition. The specification of f completes the definition of the 
C boundary condition.
C
C
C How to use the subroutine
C -------------------------
C
C The following diagram shows how the subroutine is to be used. A main
C program is required.
C
C                                   ....................................
C                                   :                                  :
C                                   :                                  :
C      ----------------------       :     --------------------------   :
C      |                    |       :     |                        |   :
C      |   MAIN PROGRAM     |------->-----|      LIBEMA            |   :
C      | (e.g. libema_t.for)|       :     |                        |   :
C      |                    |       :     --------------------------   :
C      ----------------------       :                 |                :
C                                   :                 >                :
C                                   :                 |                :
C                                   :      ------------------------    :
C          Package ---------------->:      | subordinate routines |    :
C                                   :      ------------------------    :
C                                   :                                  :
C                                   :      (this file)                 :  
C                                   :..................................:
C                                  /         |                 |
C                               |_           >                 >
C                              /             |                 |
C             ................       ................   ................  
C             :              :       :   --------   :   :  --------    : 
C             : (geom2d.for) :---<---:   | L3ALC |  :   :  | GLS |     : 
C             :              :       :   --------   :   :  --------    :  
C             :..............:       : -------------:   : -------------:  
C                                    : |subordinate|:   : |subordinate|: 
C             ................       : | routines  |:   : | routines  |:  
C             :              :       : -------------:   : -------------: 
C             : (geom3d.for) :---<---:              :   :              : 
C             :              :       : (l3alc.for)  :   : (gls.for)    :
C             :..............:       :..............:   :..............:
C                                    
C
C The contents of the main program must be linked to LIBEMA.FOR, 
C L3ALC.FOR, GLS.FOR, GEOM2D.FOR and GEOM3D.FOR.
C
C Method of solution
C ------------------
C 
C In the main program, the boundary must be described as a set of
C  panels. The panels are defined by two indices (integers) which
C  label a node or vertex on the boundary generator (the (r,z)
C  coordinate). The data structure VERTEX lists and enumerates the 
C  (r,z) coordinates of the vertices, the data structure SELV defines 
C  each panel by indicating the labels for the two edge nodes on the
C  generator and hence enumerates the panels.
C The boundary solution points (the points on the boundary at which 
C  {\phi} (SPHI) and d {\phi}/dn (SVEL) are returned) are at the centres
C  of the generator of the panels. The boundary functions {\alpha} 
C  (SALPHA), {\beta} (SBETA) and f (SF) are also defined by their values
C  at the centres of the panels.
C Normally a solution in the domain is required. By listing the (r,z)
C  coordinates of all the interior points in PINT, the subroutine
C  returns the value of {\phi} at these points in PDPHI.`
C
C
C Format and parameters of the subroutine
C ---------------------------------------
C
C The subroutine has the form
C
C      SUBROUTINE LIBEMA(MAXNV,NV,VERTEX,MAXNSE,NSE,SELV,
C     *                 MAXNPI,NPI,PINT,
C     *                 SALPHA,SBETA,SF,
C     *                 LSOL,LVALID,EGEOM,
C     *                 SPHI,SVEL,PDPHI,
C     *                 WKSPC1,WKSPC2,WKSPC3,WKSPC4,
C     *                 WKSPC5,WKSPC6,WKSPC7)

C The parameters to the subroutine
C ================================

C Geometry of the boundary S (input)
C integer MAXNV: The limit on the number of vertices of the generator
C  that defines (approximates) S. MAXNV>=3.
C integer NV: The number of vertices on S. 3<=NV<=MAXNV.
C real VERTEX(MAXNV,2): The coordinates of the vertices of the 
C  generator. VERTEX(i,1),VERTEX(i,2) are the r,z coordinates of the
C  i-th vertex.
C integer MAXNSE: The limit on the number of panels describing S.
C  MAXNSE>=3.
C integer NSE: The number of panels describing S. 3<=NSE<=MAXNSE.
C integer SELV(MAXNSE,2): The indices of the two vertices defining 
C  each panel. The generator of the i-th panel has vertices with
C  r,z coordinates (VERTEX(SELV(i,1),1),VERTEX(SELV(i,1),2)) and
C  (VERTEX(SELV(i,2),1),VERTEX(SELV(i,2),2)).

C Interior points at which the solution is to be observed (input)
C integer MAXNPI: Limit on the number of points interior to the
C  boundary. MAXNPI>=1.
C integer NPI: The number of interior points. 0<=NPI<=MAXNPI.
C real PINT(MAXNPI,2). The coordinates of the interior point.
C  PINT(i,1),PINT(i,2) are the r,z coordinates of the i-th point.

C The boundary condition ("{\alpha} {\phi} + {\beta} v = f") (input)
C real SALPHA(MAXNSE): The values of "{\alpha}" at the centres
C  of the generator of the panels.
C real SBETA(MAXNSE): The values of "{\beta}" at the centres
C  of the generator of the panels.
C real SF(MAXNSE): The values of "f" at the centres of the 
C  of the generator of the panels.

C Incident field
C real SFFPHI(MAXNSE): The incident potential at the centres of the 
C  generator of the panels
C real SFFVEL(MAXNSE): The derivative of the incident centres of the
C  generator of the panels
C real SIPHI(MAXNPI): The incident potential at the chosen interior 
C  points

C Validation and control parameters (input) 
C logical LSOL: A switch to control whether a particular solution
C  is required.
C logical LVALID: A switch to enable the choice of checking of 
C  subroutine parameters.
C real EGEOM: The maximum absolute error in the parameters that
C  describe the geometry.

C  Solution (output)
C  real SPHI(MAXNSE): The potential ("{\phi}") at the centres of the 
C   boundary panels.
C  real SVEL(MAXNSE):  ("v" or "d{\phi}/dn" where n is the outward 
C   normal to the boundary) at the centres of the boundary panels.
C  real PDPHI(MAXNPI): The potential ("{\phi}") at the interior points.

C  Working space
C  real WKSPC1(MAXNSE,MAXNSE)
C  real WKSPC2(MAXNSE,MAXNSE)
C  real WKSPC3(MAXNPI,MAXNSE)
C  real WKSPC4(MAXNPI,MAXNSE)
C  real WKSPC5(MAXNSE)
C  real WKSPC6(MAXNSE)
C  logical WKSPC7(MAXNSE)

 
C Notes on the geometric parameters
C ---------------------------------
C (1) Each of the vertices listed in VERTEX must be distinct points
C  with respect to EGEOM.
C (2) The boundary must be complete and closed. Thus 
C  SELV(i,2)=SELV(i+1,1) for i=1..NSE-1 and VERTEX(SELV(1,1),1)=0
C  and VERTEX(SELV(NSE,2),1)=0.
C (3) The indices of the nodes listed in SELV must be such that they
C  are ordered counter-clockwise around the generator of the boundary.
C (4) The generator of the largest panel must be no more than 10x 
C  the length of the generator of the smallest panel.

C Notes on the interior points 
C ----------------------------
C (1) The points in PINT should lie within the boundary, as defined
C  by the parameters VERTEX and SELV. Any point lying outside the 
C  boundary will return a corresponding value in PDPHI that is near
C  zero.

C Notes on the boundary condition
C -------------------------------
C (1) For each i=1..NSE, it must not be the case that both of SALPHA(i)
C  and SBETA(i) are zero

C External modules in external files
C ==================================
C subroutine L3ALC: Returns the individual discrete Laplace integral
C  operators. (in file L3ALC.FOR)
C subroutine GLS: Solves a general linear system of equations. 
C  (in file GLS.FOR)
C real function DIST2: Returns the distance between two 2-vectors. (in
C  file GEOM2D.FOR)

C External modules provided in the package (this file)
C ====================================================
C subroutine GL8: Returns the points and weights of the 8-point Gauss-
C  Legendre quadrature rule.
C real function FNSQRT(X): real X : Returns the square root of X.
C real function FNLOG(X): real X : Returns the natural logarithm of X.